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ABSTRACT 

Evidence for non-zero mean stellar velocities in the direction perpendicular to the 
Galactic plane has been accumulating from various recent large spectroscopic surveys. 
Previous analytical and numerical work has shown that a “breathing mode” of the 
Galactic disc, similar to what is observed in the Solar vicinity, can be the natural 
consequence of a non-axisymmetric internal perturbation of the disc. Here we provide 
a general analytical framework, in the context of perturbation theory, allowing us 
to compute the vertical bulk motions generated by a single internal perturber (bar 
or spiral pattern). In the case of the Galactic bar, we show that these analytically 
predicted bulk motions are well in line with the outcome of a numerical simulation. 
The mean vertical motions induced by the Milky Way bar are small (mean velocity 
of less than 1 km s“^) and cannot be responsible alone for the observed breathing 
mode, but they are existing. Our analytical treatment is valid close to the plane for 
all the non-axisymmetric perturbations of the disc that can be described by small- 
amplitude Fourier modes. Further work should study how the coupling of multiple 
internal perturbers and external perturbers is affecting the present analytical results. 

Key words: Galaxy: evolution - Galaxy: kinematics and dynamics - solar neighbor¬ 
hood - Galaxy: structure - galaxies: spiral 


1 INTRODUCTION 

The kinematics of stars perpendicular to the Galactic plane 
traditionally had (and still has) great importance, for in¬ 
stance as a probe of the vertical force and mass density 
in the Solar neighborhood (e.g., Kuijken fc Gilmore||1991 


Greze et al.|199^ [Siebert et al.|2003| [Bienayme et al. j2014 


Readpold l. In dynamical modeling, the natural zeroth or¬ 

der approximation is to assume a null mean vertical motion 
everywhere in the Galaxy (as it should be in a steady-state 
axisymmetric stellar system, see Binney fc Tremaine|2008 1. 
However, thanks to the spatial extension and accuracy of 
recent surveys, it has been discovered that significant non¬ 
zero mean vertical motions do exist in the extended Solar 


neighborhood (Widrow et al. 2012 Williams et al. 2013 


Garlin et al. 20131. These consist in patterns looking like 


“rarefaction-compression” waves or “breathing modes” of 
the disc. Whilst originally associated solely with excitations 
by external sources such as a passing satellite galaxy or a 
small dark matter substructure crossing the Galactic disc 


I Widrow et al.||2012 Gomez et al.|[2013 Yanny & Gardner 


2013 Feldmann & Spolyar 20151, it was then shown that 


if the density perturbation has even parity with respect to 
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the Galactic plane (i.e., is plane-symmetric), and the ver¬ 
tical velocity field has odd parity (i.e., a breathing mode). 
Then internal non-axisymmetric perturbations, such as spi¬ 
ral arms, could be sufficient to cause this (Faure et al. jMil 
hereafter F14). This was shown analytically by solving the 
linearized Euler equations in a cold fluid toy model, and 
confirmed through test-particle simulations for quasi-static 
spiral arms. Subsequently, the same effect was found in self- 
consistent simulations (Debattista 20141. Such an effect from 


internal non-axisymmetries is of course not mutually incom¬ 
patible with external perturbers playing a role in shaping the 
velocity field (e.g. Widrow et al.|[2M4 in which satellite in¬ 
teractions cause both bending and breathing modes depend¬ 
ing on the velocity of the satellite), and such perturbers are 
themselves exciting non-axisymmetric modes such as spiral 
waves I Widrow fc Bonner|2015 1. 

There is a long history of theoretical studies of the dy¬ 
namical effects of disc non-axisymmetries in two dimensions 
in the Galactic plane, dating back to the seminal works of, 
e.g., Lin & Shu (19641 and Toomre (19641 on spiral arms. 


From the observational point of view, striking kinematic fea¬ 
tures related to the bar and spiral arms were for instance 
found in the Solar neighbourhood in the form of moving 
groups, i.e. local velocity-space substructures made of stars 


of very different ages and chemical compositions (Dehnen 
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1998 [Chereul et al.|1999| |Famaey et al.|2005[|Antoja et al.| 


2008). At a less fine-grained level, it is obvious that such 


local velocity-space substructures will affect the mean mo¬ 
tions too. And indeed, it has long been well established that 
a non-zero mean radial motion is a natural response to disc 
instabilities such as spiral arms (e.g., Lin & Shu |1964| see 
also Binney &: Tremaine||200^ , independently of the exact 
nature of the perturber (quasi-static or transient). The same 
is true for the effects of the Galactic bar across the whole 
Galactic disc (see, e.g., Kuijken Tremaine|199T|. 


Such non-zero mean radial motions have recently been 
detected in the extended Solar neigbourhood with the RAVE 
survey ( Siebert et al.|20TT ), and shown to be consistent with 
the effect of spiral arms in the form of Lin-Shu density waves 
(Siebert et al. 2012). Velocity fluctuations have also been 
detected on larger scales with the APOGEE survey, and at¬ 
tributed to the effects of the Galactic bar (Bovy et al.|2015 1. 


Monari et al. (2014 hereafter M14), found that the bar can 


also explain the RAVE radial velocity gradient of |Siebert| 


et al. (2011), but did not find hints of the observed verti¬ 


cal bulk motions (Williams et al. 2013). On the contrary, 
the spiral arm model of F14 has been shown to generate a 
“breathing mode” qualitatively similar to what is observed 
for vertical motions. Generally speaking, theoretical studies 
of the dynamical response to disc non-axisymmetries in 3D 


are much less developed than in 2D (see however Patsis & 
Grosbol||1996 Cox fc G6mez|2002 ). 


In this paper, our aim is to provide a general analyt¬ 
ical framework linking the vertical bulk motions generated 
by disc non-axisymmetries to the horizontal bulk motions. 
This analytic treatment goes beyond the analytic toy-model 
of F14 for a pressureless three-dimensional fluid, which was 
only suitable for an extremely cold stellar population, and 
which required a few unrealistic assumptions (for instance 
that the vertical force of the axisymmetric potential was neg¬ 
ligible w.r.t. to the maximum vertical force of the perturber). 
Our reasoning hereafter is based on a linearization of the ze¬ 
roth order moment of the collisionless Boltzmann equation, 
i.e., the continuity equation for stellar systems. As an im¬ 
provement over the analytic treatment in F14, the present 
calculation is valid for the whole range of velocity disper¬ 
sions compatible with the epicyclic approximation, and it is 
valid close to the plane for all non-axisymmetric perturba¬ 
tions of the potential described by small-amplitude Fourier 
modes. In particular, we show that even the Galactic bar 
is expected to induce vertical bulk motions in the Galactic 
disk. We use this particular case to compare our analytic 
results with a numerical study. 


The paper is organized as follows. In Section we lin¬ 
earize the continuity equation, deriving a theoretical predic¬ 
tion for the vertical bulk motion. More specifically, we relate 
it to the horizontal motions in the case where the potential 
perturbation is a Fourier mode of small amplitude. In Sec¬ 
tion we present the outcome of a numerical test-particle 
simulation with a realistic model of the Galactic disc under 
the influence of the gravity of the Milky Way and the bar, 
and compare it with the analytical results. In Section]^ we 
conclude. 


2 VERTICAL EFFECTS OF PERTURBATIONS 
2.1 The linearized continuity equation 

Using the cylindrical coordinates {R, (j), z) and velocities 
{vr,v^,Vz) = (^R, R^, z^, and integrating the collisionless 
Boltzmann equation over velocity space, one gets the fol¬ 
lowing form of the continuity equation for stellar systems: 

dp 1 djRpuR) 1 djpu^) djpuz) ^ 

dt R dR R d4> dz ’ 

where p is the density of the system in configuration space, 
and Ui = (vi) are the average Vi velocities, all estimated at 
{R, (j), z) and at the time t. 

For an axisymmetric stationary system, symmetric with 
respect to the plane z = 0, all observables depend on {R, z) 
only, and one has 

p = po, ur^O, u^ = u^,o, Uz=0. ( 2 ) 

Let us now consider the case where the system is per¬ 
turbed by a small non-axisymmetric perturbation in the po¬ 
tential, i.e., 

■F {R, (j), z, t) = <I>o {R, z) + e<E>i [R, </>, z, t ), (3) 

where $0 is the unperturbed axisymmetric potential, e <C 1, 
and <I>i has the same order of magnitude as 4>o. 

Considering a small response to this small perturbation, 
we can write: 

Rr — eUR^ij U(j} — , Uz — eUz^i, p — poTepi. 

( 4 ) 

Plugging Eq. into Eq. 0 and dropping the terms that 
are O (e^), we obtain the well-known linearized continuity 
equation: 

dpi 1 d{RpoUR^i) po du^^i u^fidpi _ d{poUz,\) 

l)t^"R dR R dcj} Wz ■ 

( 5 ) 

From this equation, it is immediately apparent that, for a 
given background axisymmetric model, one will be able to 
relate the vertical bulk motion Uz,i to the horizontal re¬ 
sponses ur,i and u^^-i and the density wake of the pertur¬ 
bation pi. 


2.2 Solution for an exponential disc 


Disc non-axisymmetries respect the plane symmetry, hence 
on the galactic plane Uz,i{z = 0) = 0. This would not be 
true in the case where a bending mode, due to satellite in¬ 
teractions for example, is present (Xu et aL]|2015| |Widrow' 


|fc Bonnei]|2015[ ). This will be the subject of further work. 
From Uz,i{z = 0) = 0, the solution of Eq. ([^ reads 


po{R,z)uz,i{R,(f),z,t) = - R (R, 0, t) d^, (6) 


f 


where 

T! rh ^ d (RpoMj;,i) po du^^i dpi 

’ ’ ’ dt R dR R dcj) R d(j) 

( 7 ) 

We now assume that the unperturbed stellar system that we 
consider is an axisymmetric exponential disc, i.e., 


po (R, z) = Po (0,0) exp ( - 


( 8 ) 
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where /ir and are the scale length and height respectively. 
We also assume that p\ and po have the same exponential 
vertical dependence, i.e., p = pi/po is constant with 2 . With 
these assumptions Eq. § becomes 

(9) 

where 

niT) A ^ dp ^ ur,i urp ^ 9ur,i , 1 du^p 

y[K,(p,z,t)- ^ ^ OR '^R d(f> ' 

( 10 ) 


2.3 Taylor expansion close to the plane 

We now indicate only the dependence on 2 of the functions 
to simplify the notation, and Taylor expand Q in powers of 
2 up to second order close to the plane, i.e., from Eq. 


UzP (z) « - 


/;e-l^l/'^^ [g(O) + |0(o)^2]dC 


Q-\z\/hz 

Therefore, from Eq. Q and Eq. (Ill, 


( 11 ) 


Uz{z) = e sgn(z)hz 


S(0)( 


1 -e'- 


+ ^(0) (y +/i^|2| + h? - , (12) 

and, averaging Uz over 2 > 0, weighted by p, we obtain 

A/ \ ^ ^ _ o/o" dC 

A(uz} (z) = 2-^ - - 


= 2e 


fo P (0 

g(0)(hz - 2 ) + ^(O)ftz (3hz - z) 


S(0) + (6^2 + ^) 


/ - 1) 


(13) 


the difference between the mean vertical velocities of stars 
within 2 , above and below the Galactic plane, since Uz{z) = 
—Uz{—z). A value A{uz) > 0 {A{uz} < 0) implies that the 
stars of both Galactic hemispheres tend to move away from 
(towards) the Galactic plane, while for A{uz) = 0 we have 
Uz = 0, like in an axisymmetric Galaxy. 

Note that, whilst being (7° continuous, the function 


Uz(z) given by Eq. (121 is not continuously differentiable 
(C^) at 2 = 0. This is a consequence of the fact that the 
vertical density distribution exp(—| 2 |/hz) is itself not 
at 2 = 0. This technical issue disappears for a secid {z/hA 
density distribution, which is at 2 = 0. In Appendix |A| 


we compute the formulae corresponding to Eqs. (121-( 131 in 


the sech case, and show that the results are quantitatively 
similar for the same value of hz close to the plane. 


2.4 Computing the horizontal bulk motions 

At this point, we will need to estimate 0(0) and 8^0(0)/dz^. 
To do so, we consider the potential perturbation 4>i, the 
density response p, and the horizontal mean velocities as 
Fourier m-modes propagating in the disc, i.e., 

■I?! {R, (j), z, t) — Re {4>®(i?, 2 ) exp [im {4> — fist)]} . (14) 

p {R, 4>, t) = Re {p^{R) exp [im (^ — flbt)]} , (15) 


urp {R, (j), 2 , t) = Re {u%{R, 2 ) exp [im (cf) — ^^bt)j} , (16) 


Ucjtp {R, (j), z,t) = Re {u^(R, z) exp [im {(j) — flbt)]} , (17) 

where we have reintroduced the R, (j>, and t explicit depen¬ 
dencies. Dividing Eq. by Po and integrating over 2 we 
obtain the 2D continuity equation, i.e., 

dp u^fldp Urp Urp ourp 1 dU4,p 

dt^ R d4>^ R Hr ^ dR ^ R d(l> 


where 


Ui = 


Jffpiz)ui{z) d^ 
ffLPi^)dz ’ 


(19) 


and Ui is one of the urp, u^p, u%, and u^, functions 

defined above. Using the condition Eq. (181 we can derive 
p'^ as 

mU$/R + i{UyhR-UyR-dUydR) 

= - min,-U,,o/R) -■ (20) 

The perturbative regime cannot give accurate estimates 
of pUR, pu’^, and pu^^ far from the 2 = 0 plane. However, 
typically these functions must decrease quite fast as a func¬ 
tion of 2 and tend to 0 as 2 goes to infinity. Therefore, we 
can approximate 


Ui {R,<t>,t) 


jf^piR,Oui{R,<l>,U) dC 


( 21 ) 


where ^ is the height from the plane up to which we are 
computing the average and at which the integral in the per¬ 
turbative regime is converging to a constant valu^ 

We then obtain the value of u% and close the 2 = 0 
plane by linearizing Jeans equations, in the same way than 
Kuijken & Tremaine ( 1991| in the 2D case. Under the as¬ 
sumptions that the velocity dispersio ns (j\ and are re- 
lated by the epicyclic approximation (Binney & Tremaine 


20081, that the radial scale length of <l?i is larger than 


max [|itR,i|, Ur] /k, and that the mixed terms in the velocity 
dispersion or^z and are negligible, one obtains the two 
following equations for u\ and 


Ur = ^ 


^Ub — 


~Tr 


2mQ, 

R 




(22a) 


u% {R,z) = 


where 


-5$ 


a m ^mflb — 
~R 




Q.{R,z) = 


u^fi{R, 2 ) 
R 


, (22b) 


(23a) 


8fr 

K{R,zf =R^{R,z)+4n^{R,z), (23b) 

uH 


A(i?, z) = (-R, z) — z)^ , (23c) 


^ In the rest of the paper we adopt C = 1 kpc, as it is sufficient 
for the convergence of the pUi functions in the barred case of 
Section 
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B{R,z) = -n{R,z)-]^R^{R,z). (23d) 

Notice how the quantities Eq. would yield the well 

known f2, k, A, B functions ( Binney Sz Tremaine||20()8 1 on 
the z = 0 plane for a zero-asymmetric-drift cold fluid. 

The assumption that the radial scale length of "I?! is 
larger than max [|uij,i|, aij] /k is in general not valid for spi¬ 
ral arms, except for very cold populations. Nevertheless, our 
framework for relating the vertical and horizontal motions 
remains very general. In the case of spirals we should use 


the generalized equations for u% and from Lin & Shu 
||1964[), including reduction factors depending on the veloc- 


ity dispersion of the background stellar population (Binney 


& Tremaine 20081. Hereafter, we will concentrate on the 


case of the bar, where the above assumption on the radial 
variation of remains valid even for a realistic background 
stellar velocity dispersion. 


2.5 Relating the vertical and horizontal bulk 
motions 

To get the vertical bulk motions, we need G (0) and 
dG{0)/dz^ in order to evaluate A{uz){z) in Eq. (13l. To 
do so, we need and its derivatives. Hence, we further sim¬ 


plify the integrals in Eq. (211 by expanding, up to the second 


order, the functions Ui in Eq. (22aI and Eq. (22bI. These are 


even function with the respect of z, so 

Ui{R,z) « Ui{R,0) -I- 

For the specific function coming from the background 
axisymmetric model, in Appendix]^ we estimate 


(24) 


dz^ 


(R,0) 


R 

u^fi{R,0) AR 


(R), 


(25) 


where u{R) = y/d‘^^o{R,0)/dz^ is the vertical epicyclic fre¬ 
quency. 

We thus now have a fully analytical expression for 


A{uz){z) in Eq. (131, which we can apply to the case of 


the Galactic bar hereafter. Indeed, according to the above 
calculations, non-zero vertical motions should be induced 
by any non-axisymmetric perturbations described by small- 
amplitude Fourier modes, even though M14 had not found 
them in a bar simulation. This means that they are proba¬ 
bly quite small, but they should be present: here, we can ex¬ 
plicitly estimate them analytically, and then check for their 
presence in a numerical simulation similar to that of M14. 
Note however that the present analytical perturbative ap¬ 
proach breaks down close to the corotation and Lindblad 
resonances where our analytical expression for A{uz) di¬ 
verges. Note also that we consider here only the response 
to the potential perturbation, and do not relate back the 
density response to the potential perturbation. 


3 APPLICATION TO THE GALACTIC BAR 

3.1 Galactic potential and disc stellar population 

We take for the background axisymmetric potential <E>o the 
Milky Way Model I by [Binney fc Tremaine| ( |2008p . It con¬ 
sists of a spheroidal dark halo and bulge, and three disc 


components: thin, thick, and ISM disc. The disc densi¬ 
ties decrease exponentially with Galactocentric radius and 
height from the Galactic plane. The position of the Sun 
in this model is {Ro,zo) — (8 kpc,0). In this axisym¬ 
metric potential, we consider a background stellar pop¬ 
ulation described by an exponential disc of scale length 
/iR = 2 kpc and scale height hz — 0.3 kpc, with ra¬ 
dial and vertical velocity dispersions varying as ctr = 
(jR,o exp [-(R - Ro)/Rs], Uz = crz,o exp [- (R - Ro) /Ra], 
where (ctr.o, ctz^o) = (35,15) km s“^, and Rs = 5 /ir. The 
mean tangential velocity has an asymmetric drift described 
by Stromgren’s formula ( Binney fc Tremaine|2008 1 adapted 
to our case, i.e.. 


UcftfiiR, 0) = V, 


3.2 Bar potential 


2/2 

= ^ _ 


2wc 


, 7 R 

1 + F T~ 
5 hn 


(26) 


We describe the perturbation due to the bar by a 3D ex¬ 


tension of the quadrupole bar used by Weinberg (19941 and 
|Dehnen| ( |2000[ ), i.e. 'Ll as described in Eq. w with m = 2 
and the angle (f> measured from the long axis of the bar, with 


4>^(R,f) = ^( ^W(r) 


Rb 


and 


W(r) = 


-{r/RG-^ 

(r/Rb)® - 2 


for r ^ Rb, 
for r < Rb, 


(27) 


(28) 


where Ro is the Galactocentric radius of the Sun, Vo = 
'Uc(Ro) the local circular velocity, Ri, the bar length, and 
= R^ + z^. 

The amplitude e in Eq. ([^ then represents the ratio 
between the bar and axisymmetric background radial forces 
at (R, (j), z) = (Ro, 0, 0). 

Here we choose e = 0.01 (Dehnen [2000 |, Vlb = 
52.23 km s”'^ kpc“^, f2b/fI(Ro) = 1.89, ( Antoja et al.|2014 |, 
and Rb = 3.5 kpc ( Dwek et al.||1995 K 


3.3 Analytical results 

With all these inputs, we can then immediately compute 
A{uz) from Eq. (131, and in Fig. we plot the predicted 
difference between the mean vertical velocities within 300 pc 
above and below the Galactic plane, i.e. A{uz)[R,(j),z = 
0.3 kpc, t = 0). The dashed line represents the orientation 
of the long axis of the bar, and the dashed circles the position 
of the corotation and outer Lindblad resonance (OLR), i.e., 
R where B(R) = Bb, and [I1(R) — Bb] -I- k{R)/2 — 0. 

We are interested in studying the response of the back¬ 
ground axisymmetric thin disc to the bar potential, hence 
we focus on the outside of the bar region (R > Rb). We also 
limit ourselves to R < 10 kpc in this plot, since that is the 
interesting region where there are still non-negligible effects. 
In this plot the Galaxy and the bar rotate counter-clockwise. 

Fig .[^ shows that the mean vertical motions induced by 
the Milky Way bar are small, of the order of a few tenth 
of km s“^, and up to about 0.5 km s“^. More importantly, 
it shows how these motions depend on the angle from the 
long axis of the bar and the distance from the center of 
the Galaxy. The angular dependence of A{uz) is that of an 
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Figure 1. Difference between the average vertical motion of stars 
in the northern and southern hemisphere of the Milky Way A{uz) 
in our analytical model for z = 0.3 kpc, as a function of the x 
and y positions in the plane. The dashed line corresponds to the 
long axis of the bar, and the dashed circles to the position of the 
corotation and OLR. 


m = 2 mode propagating in the disk while keeping the same 
configuration w.r.t. the bar as the bar rotates. 

As noted earlier, close to the corotation and OLR, 
A{uz} diverges, because the linear theory presents poles 
there, as is evident in Eq. p0|, and Eqs. (22al-(22bl. Higher 


order expansions (and changing variables so that the new 
variables only vary slowly there) would be required near 
these resonances, where we leave our predictions as blank. 
Notice how the correspondence in the figure to the reso¬ 
nances is not exact. This is due to the fact that, since we 
consider a stellar disk rather than a cold fluid, in the ana¬ 
lytical model we use which differs from Vc for the asym¬ 
metric drift and the dependence from z, and shift the posi¬ 
tion of the resonances in the model. At the OLR there is a 
phase shift of 7r/2 in A{uz}'- along a given fixed azimuth cj), 
the value of A{uz) inside and outside the OLR has always 
opposite signs. 

At a given radius, the function A{uz) also changes sign 
between the regions ahead of the bar {(j) > 0) and behind 
the bar {(f) < 0). This is due to the fact that Uz is an odd 
function in cf), as it is the sum of odd functions, as we can find 
out by looking at Eqs. (101, (131, ( 15|17 1, ( p0| ), (22aI, and 
( |22b[ ). In particular, the compression {A{uz) < 0) is ahead 
of the bar between the corotation and the OLR, while it is 
behind the bar outside the OLR, due to the change of sign 
of A in Eq. (231. 

We note that the pattern in A{uz} is actually qualita¬ 
tively following the ur^i pattern, because they are both odd 
functions in (j>, as can be seen from Eq. (161 and Eq. (22a I 
for ur^i. This is in stark contrast with what happens in the 
case of spiral arms, for which F14 showed a clear phase shift 
between the vertical and radial bulk motions. The reason 
for this difference is that ^^{R, z) is a pure real function in 
the case of the bar, but depends on exp(i mln(R)/tanp) in 


UR(kms"h 
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Figure 2. As in Fig. [^for the horizontal motions ur (top panel) 
and Q (bottom panel) in our analytical models for 2 = 0. 


the case of spiral^ inducing a phase shift in Eq. ( |l6| and 
Eq. \22s.\ for spirals. 

For visual comparison with the A{uz) pattern of Fig. 
we plot the analytical ur and patterns for the bar 

model in Fig.[^ In particular, we see that between the coro¬ 
tation and OLR and for 0 < 0 there is a net outward motion 
of stars {ur > 0) with vertical expansion, while for (p > 0 
the radial motion is inward {ur < 0) and there is a vertical 
compression. The same behaviour has been seen in Fig. 12 
of Widrow et al. (20141, where a satellite-induced bar has 
been shown to create the same qualitative pattern. Outside 
the OLR, the pattern is reversed, both in ur and A{uz)- 

A heuristic explanation of this phenomenon is that the 
net radial flows are related to the shape of the orbits induced 


by the bar inside and outside the resonances (see Binney & 


Tremaine 20081. Outside the OLR, the orbits in the bar 


frame are aligned with the bar and retrograde, while they 
are perpendicular to the bar and retrograde between the 


2 


Where p is the pitch angle 



























6 G. Monari et al. 


corotation and OLR: this naturally explains the shape of 
the radial flow. Stars are then also adapting to the vertical 
restoring force, stronger in the inner parts than the outer 
parts of the Galaxy: when the flow is inward, the force is 
stronger and the vertical motion corresponds to a compres¬ 
sion, and when the flow is outward it corresponds to an 
expansion-rarefaction. We note that the situation is differ¬ 
ent in the case of spirals because of the restoring force from 
the spiral potential itself. 

Our analytical predictions for the vertical bulk motions 
induced by the bar can now be compared to a simulation 
akin to M14 in the next subsection. 


3.4 Numerical results 


3.4-1 Initial conditions and integration time 

A simple and yet very effective way to study numerically 
the response of a stellar disc to an external perturbation is 
through test-particle simulations. These consist of the nu¬ 
merical integration of the equations of motion of massless 
particles (representing the stars in the disc), accelerated by a 
gravitational field (representing the potential of the galaxy) 
which is uninfluenced by the particles themselves. 

We generate 5 x 10^ initial positions and velocities 
for our simulations, as discrete realizations of the Shu- 
Schwarzschild phase-space distribution function 


Bienayme & Sechaud 1997 Binney fc Tremaine] 20081 de¬ 


scribed in Section 2.2 of F14, with the same parameters 
adopted in our analytical calculations. The surface den¬ 
sity of the initial conditions obtained in this way is ap¬ 
proximately distributed in space as an exponential disc of 
scale length Iir = 2 kpc. The radial and vertical velocity 
dispersion on the disc plane vary approximately as ctr « 
(jR,oexp [-(i? - i?o)/Rs], cTz « g-^.o exp [- {R - Rp) /77s], 


where (ur.o, ctz.o) = (35,15) km s and Rb — 5/ir (Bi- 
enayme fc Sechaud|jl997[ F14). 


We then integrate forward our initial conditions for a 
total time T — 9 Gyr, from ti — —3 Gyr to te = 6 Gyr. 
For t < 0 the bar is absent and the axisymmetric gravita¬ 
tional field is only given by Model I of [Binney fc Tremaine| 
120081. This period of time allows the initial conditions 
to become well mixed in the background potentiaj^ After 
the initial 3 Gyr, we obtain a stable particle configuration, 
with velocity dispersions at {R,z) = {Ro,0)\ {(JR,cr^,az) rs 
( 37,27,13) km s“^. The vertical restoring force determines 
the 2 density profile, which, at i? = 7?o is nicely described 
by an exponential or sech^ profile (see Appendix]^, both 
with scale height ~ 0.3 kpc close to the plane. 

The bar is introduced smoothly in the simulation at 
t = 0 with the potential defined in Eq. Eq. ( [T4| | and 
Eq. (27), its amplitude reaching e = 0.01 only at t = 3 Gyr. 


For 0 ^ t < 3 Gyr the amplitude grows with time by a factor 


^ The Shu-Schwarzschild distribution function is built on approx¬ 
imate integrals of motion for an axisymmetric galaxy: in particu¬ 
lar, the vertical energy is a good approximate integral only very 
close to the plane. Because of Jeans’ theorem ( [Binney Tremaine] 
|2008| |, the worse the approximation of the integrals, the faster the 
distribution function will evolve in time. 


( Dehnen|2000 1 








The amplitude then stays constant at e = 0.01 for another 
3 Gyr. The response becomes almost stable in the rotating 
frame of the baB and the end of the simulation can be used 
to compare with our analytical predictions. 


3 . 4.3 Comparison with the analytical model 

In Fig. 1^ we present the mean vertical kinematics as a 
function of the cartesian (x, y) position in the Galaxy, at 
the end of the simulation. We plot the difference A{vz} = 
Gz)p — {vz)n, where {vz)p i{vz)n) is the average Vz of par¬ 
ticles found at 2 > 0 (2 < 0) at the end of the simulation. 
We present it both for particles at \z\ < 0.3 kpc (top panel) 
and at every 2 (bottom panel). The averages are computed 
inside bins of 0.25 kpc x 0.25 kpc, also applying a Gaussian 
smoothing on a scale 1 kpc. 

Fig-i shows a very good agreement with Fig. B with 
clear trends in the vertical kinematics, depending on the an¬ 
gle from the long axis of the bar ij) and on the distance from 
the resonances in the same way as on Fig. B The quan¬ 
tity A{vz} appears always to be periodic with (j), fluctuating 
twice for (j) € [0, 27r], so that A(vz) {<j>) ~ A{vz) {4> + t^) and 
A{vz) {<j>) ~ —A{vz} {(j) + 7^/2). Notice how the phase of the 
periodic oscillations changes inside and outside the OLR: 
the maxima of A{vz) are at ((1 ~ 7r/4 ± 7r/2 inside the reso¬ 
nance and at (() ~ —7r/4 ± 7r/2 outside, as predicted by our 
analytical model. 

The agreement between the analytical and numerical 
model is not only qualitative, but also quantitative. In 
Fig. B we show the predictions of the analytical model 
for A{uz){R,(f>,z) = A{uz){5 kpc, ((>,0.3 kpc) (blue line) 
and A{uz){R, cf>, z) = A{uz){8.5 kpc, <(), 0.3 kpc) (red line), 
where (j> is measured from the long axis of the bar. In the 
same plot we show the 90% confidence bands obtained fit¬ 
ting the model asin(2<() -|- 6) to A(vz) for particles in the 
simulation with |i? — 5 kpc| < 0.2 kpc (blue band), and 
\R — 8.5 kpc| < 0.2 kpc (red band). The concordance be¬ 
tween analytical model and simulations shows how our as¬ 
sumptions in the analytical model are convincing. Notice 
how a small phase shift is present between the analytical 
predictions and the simulation. This is due to the fact that 
the kinematics in the simulation did not reach a complete 
stability w.r.t. the bar disturbance, as previously mentioned. 


4 DISCUSSION AND CONCLUSIONS 

In this paper we demonstrated analytically how the mean 
vertical motions of stars in disc galaxies are affected by 
small, non-axisymmetric perturbations of the potential in 
the form of Fourier modes (e.g., bar or spiral arms). This 
was previously shown numerically for spiral arms in F14, 


^ Some transient effects are still present in the kinematics at the 
end of the simulation, resulting in minor asymmetries in the maps 
of the mean motions. To reach complete stability one should inte¬ 
grate for many more dynamical times, which would be unphysical 
(see, e.g., Miihlbauer fc Dehnen|2003^. 
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Figure 3. Difference between the average vertical motion of the 
northern and southern hemisphere of the simulated Milky Way 
A{v:z), as a function of the x and y positions. The dashed line 
corresponds to the long axis of the bar, and the dashed circles to 
the position of the corotation and OLR. Top panel: particles at 
\z\ < 0.3 kpc only, to be directly compared with the analytical 
model. Bottom panel: all particles. 


and supported by an analytical toy-model for a pressureless 
fluid. The present analytical treatment goes beyond this toy- 
model, and is valid for the response of realistic disc popu¬ 
lations to any non-axisymmetric perturbation described by 
small amplitude Fourier modes. 

Our general analytical solution shows that, depending 
on the distance from the Lindblad resonances and azimuth in 
the frame of the perturber, the mean vertical motion points 
towards or away from the galactic plane, with a bulk velocity 
depending on the position projected in the galactic plane 
and on the distance from the galactic plane. 

Although M14 did not find hints of significant vertical 
bulk motions in a Galactic bar simulation, our present an¬ 
alytical treatment indicates that non-zero vertical motions 
should be induced by the bar too. We thus explicitly esti¬ 
mate them analytically, and then check for their presence in 



0(rad) 

Figure 4. Predictions of the analytical model for 
A{u;^){R, <j), z) = A ( u 2)(5 kpc, (/), 0.3 kpc) (blue line) and 

A{uz)(R, <f>, z) = A{u2:}(8.5 kpc, 0, 0.3 kpc) (red line), where <f> 
is measured from the long axis of the bar. The 90% confidence 
bands are obtained fitting A(uz} for particles in the simulation 
with |ij—5 kpc| < 0.2 kpc (blue band) and |R—8.5 kpc| < 0.2 kpc 
(red band) with the model asin(20 + b). 


a numerical simulation similar to that of M14. These ver¬ 
tical mean motions induced by the Galactic bar are indeed 
modest in magnitude, especially if compared with those ob- 


served recently in the Solar neighborhood i 

Widrow et al. 

2012|Williams et al.||2013 Carlin et al.||2013 

, but they are 


well existing, as we have shown here for the very first time. 
And the results of our test-particle simulation are well in 
line with our analytical prediction away from the main reso¬ 
nances. However, it is clear that perturbations exerting ver¬ 
tical forces that change more rapidly than that of the bar 
with the distance from the galactic plane (e.g., the spiral 
arms) do create much more significant mean vertical veloc¬ 
ities. 

Our analytical treatment being valid close to the plane 
for all the non-axisymmetric perturbations of the disc that 
can be described by small-amplitude Fourier modes, it will 
be extremely useful for interpreting the outcome of simula¬ 
tions including any such perturbation, and to estimate how 
much of the observed breathing modes can be explained by 
such non-axisymmetries alone. Further work should study 
how the coupling of multiple internal perturbers (bar-1- mul¬ 
tiple spirals) and external perturbers (satellite interactions 
themselves exciting spiral waves, but also bending modes) 
is affecting the present analytical results. 
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APPENDIX A: RESULTS FOR A SECH^ 
VERTICAL DENSITY DISTRIBUTION 


We report here the analytical results for and A{uz} ob¬ 
tained with a different unperturbed vertical density distri¬ 
bution, namely 



sech^ 



(Al) 


With this continuously differentiable choice of density, and 
following the same steps as in Section]^ Eqs. (12l-(13l be¬ 
come 


^{z) = —ehz cosh^ 


Q(0) tanh 


-s<»> 


hw 


24 


hz 

-2z/hz 


— hz In 1 J- e' 


-2z/h: 




^ - 1 u 


(A2) 


and, 

A{uz) = —2e coth ( -A ] -j —hzC/(0) In 


cosh ( — 
hz 


-s<»> 


4 


9h® 

—C(3)- —LU (- 

16 sw 2 

..,2 1 2 
TT 

24 


Lig ^-e 

(-efc) 


Z+^Ll2 (- 


i2(-e 

2 + Y log (e"ft -f l) 2^ -b y j. 


(A3) 


where the “polylogarithm” function Lis ( 2 ) is defined by the 
power series 


00 ^ 

Li42) = E|7’ (A4) 

k=l 

and the Riemann (^(s) function by 

00 

C(s) = Efc''- (A5) 

fc=l 


corresponding to Fig. and Fig. 1^ corresponding to 
Fig. El These figures show how in the bar case the rela¬ 
tive difference between the predicted bulk motions for the 
case of a sech^( 2 /hz), and exp(—|2|/hz) density distribu¬ 
tions are tiny (of the order of 5 x 10“^ km s“^ at most). 
We note that the initial axisymmetric stellar population in 
our test particle simulation is well fit close to the plane 
by sech^ and exponential distributions with the same scale 
height hz = 0.3 kpc. However, we note that the sech^ pro¬ 
file is actually a better fit, as expected from the form of the 
Shu-Schwarzschild distribution function used for the initial 
conditions. 


Using these new formulae Eqs. (|A2l)-( A31, we obtain Fig. Al 


APPENDIX B: ESTIMATING u^,o NEAR THE 
PLANE 


tern, and neglecting mixed terms can be rewritten as (Binney 
& Tremaine|200'8 1 


9 (Poo-|,o) 
dR 


+ Po 






— U^,0 


R 


d^o 

dR 


= 0. (Bl) 


Differentiating twice with the respect of 2 , assuming an ex¬ 
ponential disc, and that the velocity dispersion near the 
plane is approximately constant with z, we obtain 


mIo {R, z) + M0,o {R, z) ^ Q^2° (-^’ 


R 

2 dRdz"^ 


{R,z), 

(B2) 


po {R, z) = Po (0, 0) exp 
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Figure Al. As in Fig.^ but for a density distribution depending 
on 2 as sec\d{z/hz). 



0(rad) 

Figure A2. Comparison between the analytical A{uz) in the 
case of a density distribution depending on z as eicp{—\z\/hz) 
(colored lines) and sec\d(z/hz) (black dashed lines). The poten¬ 
tial parameters and scale height are as in Fig. The blue line 
correspond to H = 5 kpc, the red line to R = 8.5 kpc. 


which, estimated at 2 = 0, is 


and 


dz^ 


(i?,0) 


R dv"^ 
2u4,fi{R,0) dR 


(R), 


v\R) 


dz^ 


{R,0), 


is the square of the vertical frequency. 


(B3) 


(B4) 











